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ABSTRACT 

We address the problem of identifying remnants of satellite galaxies in the halo of our galaxy 
with Gaia data. The Gaia astrometric mission offers a unique opportunity to search for and 
study these remnants using full phase-space information for our Galaxy's halo. However, the 
remnants have to be extracted from a very large data set (of order 10 9 stars) in the presence 
of observational errors and against a background population of Galactic stars. We address this 
issue through numerical simulations with a view towards timely preparations for the scientific 
exploitation of the Gaia data. 

We present a Monte Carlo simulation of the Gaia catalogue with a realistic number of 
entries. We use a model of the galaxy that includes separate light distributions and kinematics 
for the bulge, disc and stellar halo components. For practical reasons we exclude the region 
within Galactic coordinates: —90° ^ £ ^ 90° and —5° ^6^5°, close to the Galactic plane 
and centre. Nevertheless, our catalogue contains 3.5 x 10 8 stars. No interstellar absorption 
has been modelled, as we limit our study to high Galactic latitudes. 

We perform tree code 10 6 -body simulations of satellite dwarf galaxies in orbit around 
a rigid mass model of the Galaxy. We follow the simulations for 10 10 years. The resulting 
shrinking satellite cores and tidal tails are then added to the Monte Carlo simulation of the 
Gaia catalogue. To assign photometric properties to the particles we use a Hess diagram for 
the Solar neighbourhood for Galactic particles, while for the dwarf galaxy particles we use 
isochrones from the Padova group. When combining the Milky Way and dwarf galaxy models 
we include the complication that the luminosity function of the satellite is probed at various 
depths as a function of position along the tidal tails. The combined Galaxy and satellites model 
is converted to a synthetic Gaia catalogue using a detailed model for the expected astrometric 
and radial velocity errors, depending on magnitude, colour and sky position of the stars. 

We explore the feasibility of detecting tidal streams in the halo using the energy versus 
angular momentum plane. We find that a straightforward search in this plane will be very 
challenging. The combination of the background population and the observational errors will 
make it difficult to detect tidal streams as discrete structures in the E-L z plane. In addition 
the propagation of observational errors leads to apparent caustic structures in the integrals of 
motion space that may be mistaken for physical entities. Any practical search strategy will 
have to use a combination of pre-selection of high-quality data and complementary searches 
using the photometric data that will be provided by Gaia. 

Key words: Galaxy: formation, structure - galaxies: interactions - methods: numerical 



1 INTRODUCTION 

Current cosmological models for the formation of the large scale 
structure of the Universe envisage the formation of large lumi- 
nous galaxies as a long process of merging of smaller structures 
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( Kauffman n & Whitelll993l) . Although the details that lead from 
the merging of dark haloes to the formation of luminous galax- 
ies are quite complicated and currently addressed only through the 
use of semi-analytical mod els l Kauffmann, White & Guiderdoni 
1 19931 ISomerville &Primacldll999l:IColeatalJl2000l) . it is clear 
that such a formation scenario must have left an imprint on the 
halo of luminous galaxies like ours, where the long dynamical 



2 A.G.A. Brown, H.M. Velazquez & L.A. Aguilar 



Table 1. Scientific capabilities of the Gaia mission. Listed are the survey parameters and the measurements with 
their expected accuracies. Note that ^arcsec stan ds for micro-arcsecond. T he numbers in this table are from the 
Gaia Concept and Technology Study Report (see ESA 2000; Perrvman 2002) 



Survey parameters Measurements and accuracies 



Magnitude limit 


20-21 mag 


Astrometry 


4 /^arcsec at V = 10 




Completeness 


20 mag 




1 1 /^arcsec at V = 15 




Number of objects 


26 million to V = 15 




160 /iarcsec at V = 20 






250 million to V = 18 


Photometry 


4 broad band to V = 20 






1000 million to V = 20 




1 1 medium band to V = 


20 


Observing programme 


On-board and unbiased 


Radial velocities 


1-10 km s" 1 at V = 16 


-17 



tim e scales preserve fo r a long time the remnants of old merg- 
ers <Lvnden-Belllll976l) . Although with low statistical weight due 
to the small number of stars involved, for a long time there 
has been increasing evidence of the existence o f substruc ture in 
the st el lar halo of our galaxy |E ggen|ll965{: [Rodgers & Paltogloul 
[ 1984 iRatnatunea & FreemarJ 1198.4 iDionidis & Beers! Il989t: 
lArnold & Gilmorelll992tlMaiewski. Munn & Hawlevll 19961) . This 
trend culminated with the discovery of the Sagittarius Dwarf galaxy 
bv llbata. Gilmoj^^r^in| ^^4|) and the discovery of the Canis 
Major dwarf I Newber g et alj2002tlMartin et alj20 041. These satel- 
lite galaxies are in the process of being torn apart by the tidal field 
of our Galaxy. 

The identification and study of substructure in our stellar 
halo is of paramount importance, not only because it traces the 
formation of our Galaxy, but because it allows us to probe the 
clustering of da rk matter at small spatial sca les (for a recent re- 
view see Freeman & Bland-Hawthorn 2002). This clustering is 
a crucial clue to the nature of its constituents: The nature of 
the dark matter sets the scale at which dark matter can clump 
teond&Szalavll983tlN avarro. Frenk & White! 1997tlKTvpin et alJ 
ll999UMooreet alJll9991) . A popular current candidate, the axion, 
is "cold" enough that it can clump at scales of dwarf galaxies. Al- 
though evidence for dark matter within dwa rf galaxies orbiting our 
galaxy has been put forward llrwin & Hatz idimitriou, fl995l) . this 
evidence is weakened by the unfortunate degeneracy between po- 
tent ial well depth and orbital eccentricity in the dynamical mod- 
els 1 Wilkinson et al. 2002). Identification and modelling of several 
satellite systems in the process of being torn apart by the tidal field 
of our galaxy can help to disentangle this degeneracy, as the rate of 
erosion of the satellite depends more importantly on the depth of 
its own potential well t han on the details of its internal dynamics 
(Roche Lobe argument, Mihalas & Binney 1981). 

Several authors have pointed out that the thick disc of our 
galaxy may have been pro duced by a merger even t around 1 
Gyr a g o jGilmore & Wvsdll985|: ICarnev. Latham & Lairdlll989h . 
Wvse ( 1999) has pointed out that this corresponds to a red- 
shift of z ~ 2. By probing the structure that exists within 
the halo of our galaxy, we are doing 'in situ' cosmological 
studi es of events that took place long ago. Recent ob serva- 
tions IStockton. Canalizo & Maihara 2004 ; Tecza et al. 2004J) have 
shown evidence for the existence of luminous galaxies at large 
redshifts (z ~ 2.5) that seem to be old, massive and have large 
metallicity. Their existence in large numbers, if confirmed, poses a 
challenge to the standard cosmological scenario of bottom-up for- 
mation. It is clear that it is of paramount importance to study the 
stratigraphy of the halo of our galaxy to try to discern what the 
assemblage process was by which our galaxy arose. 

Great progress has recently been made in the study of sub- 
structure in the halo of the Milky Way through dedicated surveys 



of ha l o giants and other tracer populations (e .g., iMorrison et alJ 
1200(1 iMaiewski et ai]|2000l : IPalma et aljl2003h . In ad dition large 
scale photometri c surve ys, such as SDSS and 2MASS lYork et aU 
teOOdlCutri et alj|2003t) which cover the entire sky or large parts 
thereo f, have been exploited to explore t he halo substructur e (e.g ., 
iNewberg et al1l2002t [Martin et alJboollRocha-Pinto et alJ 2004). 
However all these surveys provide at most one component of the 
space motion for a very small fraction of the stars surveyed and 
often they only probe a small fraction of the halo. To definitively 
unravel the structure and formation history of the Galactic halo (and 
the Galaxy as a whole) detailed phase-space information is required 
of a sig nificant fracti on of the enti re volume of the Galactic halo 
jFreeman & Bland-Hawthornl2002t . This requires a very high pre- 
cision astrometric survey done over the whole sky, complemented 
with radial velocities and photometry. 

In 2000 October the Gaia satellite mission was selected as 
the European Space Agency's sixth Cornerstone mission to be 
launched no later than 2012. The main scientific aim of the Gaia 
mission is to map the structure of the Galaxy and unravel its for- 
mation history by providing a stereoscopic census of 1 billion stars 
throughout the Galaxy. This will be achieved by measuring accu- 
rate positions, parallaxes and proper motions for all stars to a lim- 
iting magnitude of V = 20. To complement the astrometric infor- 
mation radial velocities will be collected for all stars brighter than 
V ~ 17. In addition photometry will be measured for all stars; four 
broad bands will be employed and about 11 intermediate bands, 
which will provide detailed astrophysical information for each de- 
tected object. Gaia will employ an on-board detection scheme al- 
lowing it to measure every detected object brighter than V = 20. 
This ensures an unbiased and complete survey of the Galaxy. The 
scientific capabilit ies of Gaia are summarized i n Tabled For de - 
tails see lPerrvmar]<2002l) . Ferrvman et ai]<200ll) . and lESAli200oT) . 

These measurement accuracies can be put into perspective 
by using Galaxy models to calculate the corresponding distance 
and tangential velocity accuracies (see Perrvman 2002): 21 million 
stars will have distances with a precision better than 1 per cent, 116 
million better than 5 per cent and 220 million better than 10 per 
cent; 84 million stars will have tangential velocities with a precision 
better than 1 km s ~ 1 , 2 1 million better than 3 km s _ 1 , 300 million 
better than 5 km s , and 440 million better than 10 km s _1 . 

In view of this kind of stereoscopic data becoming avail- 
able for the Galaxy in the future, several authors have made 
simulations of the process of accretion and tidal disruption 
of satellite galaxies in the halo of our galaxy, with the 
goal of identifying the ch aracterist ics of the resulting debris 
(e.g.lOh. Lin & Aarsetr]|l995Ujohriston Hernquist & Boltdll996t 
iMeza et all2005l) . lHelmi & WhiteNl999l) . and lHelmi & de Zeeuwl 
(2000), in particular, have studied the problem of identifying the 
resulting substructure in the stellar halo with information from a 
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new generation of astrometric space missions. They conclude that 
a promising technique is to search for this substructure using con- 
served dynamical quantities: energy and angular momentum. They 
performed simulations of the destruction and smearing of satellite 
galaxies in the halo of our galaxy and plotted the remnants in the 
(E, L, L z ) space. They applied a simple model of the astrometric 
errors expected in the FAME, DIVA and Gaia missions, that de- 
pends on the apparent magnitude alone. The errors for the Gaia 
mission will depend also on the Ecliptic coordinates and the colour 
of stars, and this can introduce additional correlations between the 
astrometric variables. In addition [Helmi & de Zeeuwl feOOd) only 
considered relatively bright and nearby (V < 15, d < 6 kpc) 
stars when constructing their E-L-L z diagrams. However, when 
searching for substructure throughout the whole volume occupied 
by the Galaxy's halo, the majority of stars will be much fainter (and 
much further away) and have correspondingly worse astrometric er- 
rors (see Table0 which will result in much more smearing in the 
(E,L,L Z ) space. Finally, previous investigations have not taken 
into account the problem of contamination by the Galactic back- 
ground. It is expected that, in any given field, the signature of an 
accretion event will not dominate the corresponding star counts and 
a method should be devised to identify the proverbial 'needle in the 
haystack' . 

The goals of this work are to: 

(1) understand, through simulations, how one can detect and 
trace the debris of disrupted satellites in the halo of our galaxy using 
direct measurements of phase-space, as will be provided by Gaia. 

(2) include a realistic model of the Galactic background popula- 
tion against which the remnants have to be detected. This includes 
a more elaborate model of the observational errors expected for 
Gaia. 

(3) gain practical experience with analysing and visualizing the 
enormous volume of multi-dimensional information that will be 
present in the Gaia data base. 

For the Galaxy model we have been careful to use Monte 
Carlo simulations that result in the expected number of stars in the 
Gaia catalogue to avoid problems owing to under sampling. Our 
Galactic model includes spatial, photometric and kinematic infor- 
mation, as all these factors affect the final astrometric precision. For 
the simulated satellites, we tackle the problem of probing the satel- 
lite luminosity function at varying depths as a function of distance 
to the observer, as we move along the tidal streamer. We present a 
method to avoid loosing too many iV-body particles from the result- 
ing simulated survey, while keeping a realistic model of the depth 
of the observations. At the present time we have not incorporated 
the effects of dynamical friction, and this remains the main fac- 
tor not yet t aken into account. Simulations like those presented by 
iMeza et suggest that this may be an important effect that 

further smears the expected signal in the E versus L z plane. We 
plan to incorporate this effect in the near future, however, our in- 
corporation of the Galactic background, the more realistic astrom- 
etry errors and the addressing of the sampling details, introduces 
enough new elements that publishing these results is warranted. 

In Section|2|we present our photometric and kinematic model 
for the luminous components of the Galaxy and discuss the prob- 
lems that have been tackled to ensure a realistic modelling of the 
magnitude limited probing of Gaia. Section [3] presents the details 
of the A-body simulations that we performed to study the process 
of tidal dissolution of the satellites orbiting in the potential of the 
Milky Way galaxy. Section[4]presents the simulated Gaia survey in 
detail, describing the model used for the observational errors and 



the way the A-body information has been incorporated into the 
simulated survey. In Section [5] we discuss the appearance of the 
satellites in the E versus L z plane, emphasizing the effects of in- 
cluding a background population and a more detailed error model. 
Finally, Section|6|presents our conclusions and directions for future 
work. 



2 MONTE CARLO MODEL OF THE GALAXY 

The Galaxy simulations aim at providing a synthetic Gaia cata- 
logue consisting of all stars on the sky brighter than V = 20. This 
catalogue will consist of a very large number of stars (~ 10 9 ) com- 
prising a mix of various populations. The main Galactic compo- 
nents will be stars in the bulge, disc and halo. Superposed will be 
the satellites presently orbiting the Galaxy (such as the Magellanic 
Clouds) as well as debris streams of satellites that have been ac- 
creted and have not yet completely mixed with the population of 
Galactic stars that was in place before the accretion events. 

We have simulated in detail the disruption of satellites on var- 
ious orbits around the Galactic centre, assuming a smooth Galac- 
tic potential consisting of a disc, bulge and halo component (see 
Section|3j. This smooth component corresponds to the above men- 
tioned population of Galactic stars that was in place before the ac- 
cretion event started. It is against this (vastly larger) 'background' 
population of stars in phase-space that one has find the tidal de- 
bris. One of our main concerns here is a realistic simulation of this 
background in the Gaia catalogue. 

We can approach this problem in two ways, (a) We can assume 
that the populations of stars constituting the bulk of our galaxy is 
already thoroughly mixed, both spatially and kinematically, so that 
they form spatially smooth components describable with simple 
kinematics (rotation plus velocity dispersions). Here we implicitly 
assume that the bulk of our galaxy was formed very early on in 
some form of monolithic collapse or that the traces of any mergers 
resulting in the formation of the Galaxy have been wiped out in the 
phase mixing process, (b) The hierarchical paradigm for galaxy for- 
mation assumes that the first structures to form in the Universe are 
small galaxies which through many mergers give rise to the large 
galaxies such as our own. Depending on the merger history of the 
Galaxy it could thus well be that the halo is not smooth but consists 
of a superposition of many merged small galaxies which will lead 
to a lumpy phase-space structure. In this case the bulk of the Galaxy 
would consist of smooth bulge and disc populations but a halo in 
which the stars are clumped both spatially and kinematically. 

These two options for modelling the background populations 
will translate in to different distributions of the stars in phase-space. 
We want to stress here that our main concern is not to build an ac- 
curate model of the Galaxy but to use the model we have to obtain 
a distribution of stars in phase-space that captures the essential fea- 
tures of any population that has been in place for some time before 
the accretion events that we wish to find. These features are (1) 
that these stars occupy the entire volume of the Galaxy according 
to a density distribution that is much smoother than that of the de- 
bris of recent accretion events; and (2) that these stars have mixed 
to a degree where their kinematics is much simpler (rotation plus 
dispersion) than that of the accretion debris. 

Thus, although model (b) seems to be supported by recent ob- 
servational evidence we chose model (a) for simplicity. Certainly 
in the outer reaches of the Galactic halo this may lead to unreal- 
istically smooth distributions of stars in phase-space. However, an 
analysis by New berg et alj 12002. based on SDSS data) of the spa- 
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Table 2. The three spatial components of the Monte Carlo model of the 
luminosity distribution of the Galaxy. Their analytical forms (except for the 
normalizations) are given as well as the length scales involved (in kpc). 
The variables r and (R, z) indicate spherical and cylindrical coordinates 
respectively, centred on the Galactic centre. 

Bulge (r|+r 2 )~ 5 / 2 r B = 0.38 

Disc exp (-( " B + M)) R Q =8.5 R u = 3.5 

z D = 0.2 

Halo (r 7 / 2 + (r') 7/2 )~ 1 r H = 1 <?h = 0.8 

r> = (i? 2 + (^/ te ) 2 )V2 



tial distribution of stars near the turn-off in the Hertzsprung-Russell 
diagram clearly shows the presence of a smooth component of the 
halo. Moreover, in the inner parts of the Galaxy (near the Sun's po- 
sition) the approximation of our Galaxy as consisting of a smooth 
set of stellar populations will be valid even if the Galactic halo con- 
sists of superpositions of streams of many accret ed satellites. This 
was argued recently bv lHelmi. White & Springel l2002t) on the ba- 
sis of the analysis of high-resolution cosmological simulations of 
the formation of galactic haloes. They showed that in these simu- 
lations the phase-space distribution of dark matter particles in the 
Solar vicinity is very smooth spatially, with the velocity ellipsoids 
deviating only slightly from Gaussian. 

Finally we note that one could use these same cosmological 
simulations to build a more realistic Galaxy model, although one is 
then faced with the problem of translating the phase-space param- 
eters of ~ 10 7 dark matter particles into those of ~ 10 9 stars. 

We use two different Galaxy models in this paper. The Monte 
Carlo model generates a distribution of stars (i.e., luminous objects) 
for the three components of the Galaxy, and these stars make up the 
bulk of the entries in our simulated Gaia catalogue. We will refer 
to this model as the 'light' model of the Galaxy. For the iV-body 
simulations of the satellites we need a Galactic model that includes 
both the luminous and dark matter. This model is represented by a 
gravitational potential only (see Section |5J and we will refer to it 
as the 'mass' model. We now discuss the Galactic light model and 
how we implemented its Monte Carlo realization in practice. 



2.1 The Galactic luminosity distribution model 

The model for the light distribution of the Galaxy consists of three 
spatial components; a bulge with a Plummer law density distri- 
bution, a double exponential disc and a flattened halo component 
for which the density falls of as r -3 ' 5 at large distances from the 
Galactic centre. The analytical forms together with the scalelengths 
of these components are listed in Table|2| The component luminosi- 
ties are all normalized to their luminosity density at the position of 
the Sun. We assume a total luminosity density at the position of 
the Sun of 6.7 x 10 -2 L© pc~ 3 and bulge-to-disc and halo-to-disc 
luminosity density ratios of 5.5 x 10~ 5 and 1.25 x 10™ 3 , respec- 
tively. This corresponds to a bulge-to-disc luminosity ratio of 0.20 
and a halo-to-disc luminosity ratio of 0.1667. The corresponding 
total luminosity is 3.2 x 10 10 L©. 

The stars in this model are all assigned a spectral type and lu- 
minosity class according to the Hess-diagram listed in Table 4—7 of 
Mihalas & Binne^ <198ll) . The Hess-dia gram provides the relative 
numbers of stars in bins of absolute magnitude (My) and spec- 
tral type. These numbers integrated over spectral type provide the 



Table 3. The kinematics of the three components of the model Galaxy. The 
rotation velocity and velocity ellipsoid(s) of each component are listed in 
units of km s —1 . The velocity ellipsoid of the bulge is isotropic. For the 
disc and halo the velocity ellipsoids are given in Galactocentric cylindrical 
coordinates as function of stellar spectral type. 



Component 


Vrot 


Velocity ellipsoid 






Bulge 





a = 110 








Disc 


220 


Sp. Type 


a"/? 




c z 






O 


10 


9 


6 






B 


10 


9 


6 






A 


20 


9 


9 






F 


27 


17 


17 






G 


32 


17 


15 






K 


35 


20 


16 






M 


31 


23 


16 


Halo 


35 


(tr = 135, 




= 105, 


<r z = 90 



luminosity function. We consider the Hess-diagram fixed for all 
Galactic components and all positions throughout the Galaxy. 

Finally, the kinematics are modelled assuming that each 
Galactic component rotates with a constant velocity dispersion. For 
the bulge an isotropic velocity dispersion is assumed. For the disc 
a different velocity ellipsoid is assumed for each spectral type. For 
the halo the same velocity ellipsoid is assumed regardless of spec- 
tral type. In all cases, the velocity ellipsoid is assumed to be fixed 
in cylindrical coordinates. The velocity dispersions and rotation ve- 
locities are listed in Table [5] Th e values for the bulge h ave been 
taken from section 10.2.4 in Binnev & Merrifieldl il998) (except 
for the rotation which we have arbitr arily set to zero). Th e values 
for the disc are from table 7-1 in Mihalas & BinnevHl98lh (ignor- 
ing the vertex deviati ons of the ve locity ellipsoids) and for the halo 
table 10.6 in Binn ev & Merrifield (1998) was used. 

Although this is a highly simplified model of the Galaxy it is 
good enough for providing the desired 'background' distribution in 
phase-space. 



2.2 Practical implementation of the Monte Carlo model 

A straightforward realization of our Galaxy model consists of inde- 
pendently generating for each star a random space position, abso- 
lute magnitude and population type, followed by the corresponding 
kinematics. However, the Gaia survey will be magnitude limited 
and in that case this simple Monte Carlo procedure is potentially 
very wasteful. For each star one has to check whether it will be 
bright enough in apparent magnitude to be included in the survey. 
If the magnitude-limited volume is small compared to the volume 
occupied by the Galaxy, a large fraction of stars will not appear in 
the final survey. This is a particularly severe problem for intrinsi- 
cally faint stars, which are the most abundant. 

The strategy we use minimises wasted effort by generating 
stars only within the magnitude limited sphere. Let us assume a 
Galaxy whose density (in number of stars per unit volume) is given 
by pG (r) and whose luminosity function (fraction of stars with ab- 
solute magnitudes between My and Aiy + 1) is given by 4>g(Mv), 
regardless of position in the Galaxy. An observer at the position 
?"obs within the Galaxy conducts a magnitude limited survey with a 
cut-off at mii m . This means that for stars with absolute magnitudes 
between My and Mv + 1 the number of stars in the survey will 
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be: 



rK is {M v ) = / p G (r)0 G (Mv)dr 



(1) 



= 4>g(M v ) / p G (r)dr, 

where Q.v is the magnitude limited sphere centred on the observer, 
whose radius is given by the maximum distance at which the star is 
still visible: 



0.2(m lim -M v ) + l 



(2) 



The total number of stars can be found by integrating over all ab- 
solute magnitudes: 



JV& = / n° s (Mv)dM v 



</>g(M v 



p G (r)dr 



dM v 



(3) 



<j>G(Mv)t G (Mv)dMv , 



where M v and M v are the faint and bright limits of the luminosity 
function and Qg(Mv) is the integral within square brackets which 
represents the number of stars in the Galaxy within the magnitude 
limited volume Qv- Taking the product of £g(Mv) with (/>g(Mv) 
brings down (pG(Mv) to the actual number of stars of magnitude 
between My and Mv + 1 seen in the survey. 

If we generate a Monte Carlo simulation of the magnitude 
limited survey with N^c 'stars' such that the luminosity func- 
tion is randomly sampled and the density of the galaxy is fully 
and randomly sampled within the Qv corresponding to each Mv, 
the Monte Carlo simulation will represent a 'fair' sample of what 
would be obtained in the magnitude limited survey. 

One way of ensuring this proper random sampling is as fol- 
lows. We define a weighted luminosity as: 



^c(Mv) = 4>g{M v )C,g{Mv) , 



(4) 



and draw a random Mv from iPg(Mv)- Subsequently we draw a 
random position from within the corresponding magnitude limited 
sphere Qv, using pG as a probability function within this volume 1 . 
We take an upper limit to the radius of the magnitude limited vol- 
ume of 100 kpc. In Fig.Q we show both the intrinsic and weighted 
luminosity functions <pG and %f)G for our galaxy model. The in- 
trinsic luminosity func tion was obtained by sum ming over spectral 
type in Table 4-7 of iMihalas & Binned <198ll) and the weighted 
luminosity function was calculated for — 5 ^ Mv ^ +16. 

The Monte Carlo realization of our galaxy model thus pro- 
ceeds as follows. Given that the luminosity function and spectral 
type distribution is assumed the same for all Galactic components 
and positions within the Galaxy, we start by drawing a random 
value for Mv and spectral type from the Hess diagram, weighted in 
luminosity according to the procedure outlined above. The values 
for Mv are only tabulated as integers which will produce disconti- 
nuities in the the generated stellar positions. To avoid this we add a 
random fraction between and 1 to each generated value of My. 

Once the absolute magnitud e of a star is chos en, the Von Neu- 
mann rejection technique (e.g. IPress et allll992L Section 7.3) is 
used to draw a random 3D position within the corresponding mag- 
nitude limited sphere. The sum of the densities of each Galactic 

1 this means a renormalization such that p G = 1 



ui_ 4 




Figure 1. Intrinsic and weighted luminosity functions <f>G (left vertical 
scale) and i/>g (right vertical scale) for our galaxy model. The circles indi- 
cate 4>G ar >d the triangles ipQ , plotted as the actual number of stars expected 
in the survey. 



component, with the appropriate relative amplitudes, is used for 
this task. Once the position is obtained, the Galactic component 
to which the star belongs is assigned using the relative amplitude 
of each component at that position. Finally, a velocity vector is 
randomly drawn from the kinematical model appropriate for the 
Galactic component to which the star has been assigned. 

We use the Von Neumann rejection technique to generate the 
positions because the total density distribution is a rather compli- 
cated function of the cylindrical coordinates R and z when these 
are centred on the observer. Although simple to implement, this 
technique turns out be very inefficient for our case. We want to 
generate stellar positions within the magnitude limited sphere Qv 
centred on the position of the Sun. To apply the rejection technique 
we need a functio n which is every where within Qv larger than pa 
(see section 7.3 of Pr ess et alll992l) . The simplest function to use is 
a constant equal to the maximum density within this sphere, how- 
ever its value cannot be obtained in a simple way. So one is then 
forced to use the density at the Galactic centre as the upper bound 
to pa and then first generate a position for the star in Qv before 
applying the Von Neumann technique to decide on whether or not 
to retain the star. For positions far away from the Galactic centre 
(either faint stars near the Sun or brighter stars far out in the halo) 
this becomes very inefficient owing to the large density contrast. 

To prepare for the volume of data that Gaia will deliver and in 
order to have a realistic number of Galactic background stars with 
respect to the satellite stars, we want to achieve a one-to-one sam- 
pling of the Gaia survey in our Monte Carlo model of the Galaxy. 
This implies that we need to generate the positions of ~ 1.5 x 10 9 
stars. This is a task that would take a prohibitive amount of time, 
even using a fast cluster of Pentium CPUs as we did (see Section|3j. 

In practice the tracing of satellites on the sky in the parts of 
the Galactic plane towards the inner Galaxy will be very hard ow- 
ing mainly to the large extinction in those directions. Therefore we 
decided not to simulate the part of the sky with Galactic longitude 
between —90° and +90°, and with Galactic latitude between —5° 
and +5°. In our Galactic model ~ 80 per cent of the stars lie in 
this region of the sky as seen from the Sun for a survey limited at 
V = 20. Hence we now have to simulate 'only' ~ 3 x 10 8 stars. 
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Table 5. Self-consistent satellite models (10" particles). 



Figure 2. Rotation curve for our Milky-Way mass model used in the TV- 
body simulations. 



Table 4. Parameters for our Milky Way mass model used in the TV-body 
simulations. 



Disc 




Bulge 


Halo 


M d = 


5.6 x 10 10 Mq 


M b = 1.4 x 10 10 Mq 


Vh = 186 kms — 1 


Rd = 


3.5 kpc 


a = 630 pc 


Rc = 12 kpc 


P~ 1 = 


= 700 pc 




<?* h = 0.8 



Moreover, cutting out this part of the sky (and thus the Galactic 
centre) allows us to use a much smaller upper bound for the den- 
sity which turns out to be the biggest contribution to improving 
efficiency. Now the entire Monte Carlo Galaxy can be generated in 
two weeks' time on a Beowulf computer. 

The resulting simulated survey of the Galaxy is thus fully sam- 
pled and contains 3.5 x 10 s stars. For each of these we generated 
the 6 phase-space coordinates, an absolute magnitude, a popula- 
tion type (bulge, disc or halo) and a spectral type. In Section|4]we 
discuss in detail how these data were converted into a simulated 
Gala survey consisting of positions on the sky, parallaxes, proper 
motions, radial velocities and photometric information. 



3 TV-BODY SIMULATIONS 

In this section we briefly describe the numerical technique used to 
study the disruption of satellites by the Milky Way's gravitational 
potential. Given that our aim is to identify the satellite remnants, 
we have chosen a rigid potential to represent the Milky Way while 
self-consistent TV-body realizations are employed for the satellites. 
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3.1 Numerical models 

3.1.1 Galaxy mass model 

Our Galaxy mass model for the TV-body simulations is a composite 
consisting of a disc, a bulge and a halo. The axisymmetric disc 
component is represented by the following potential-density pair 
(Quinn & Goodman 1986): 



Pd(R,z) = 



Mdj 



(5) 



GMd 



dk 



-fc|z| 



kJ {kR) (3 2 
(k 2 + /3 2 ) 3 / 2 (3 2 - k 2 
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where Md is the mass of the disc, Rd is the radial scalelength and 
/3 _1 is the vertical scaleheight. 

A spherical Hernquist profile has been adopted for the bulge 
component (Hernquist 1990): 

Mb a 



Pb{r) = 



2n r(r + a) 



(6) 



where Mb is the bulge mass and a its scalelength. Another key in- 
gredient of our model is the massive dark halo which is represented 
by a logarithmic potential iBinnev & Tremainell987|) . 



<S> h = v 2 ln(R 2 c +R 2 + z 2 /ql h ) 



(7) 



where Vh, Rc and q$ h = (c/a)$ h are constants. The last param- 
eter is the flattening of the halo potential. Table |4] summarises the 
parameters of our Milky Way mass model and its rotation curve is 
shown in Figure|2| The total mass of the Milky Way model, trun- 
cated at 200 kpc, is 1.7 x 10 12 M . 

Current models of large scale structure and galaxy formation 
postulate that small objects are t he first to form and t hen aggregate 
into larger systems (e.g. iKauffmann & Whitelll993l) . In the non- 
baryonic CDM model, this hierarchical assembly of haloes leads 
to triaxial shape s with average axis ratios o f (c/a) Ph = 0.5 and 
(b/a) Ph — 0.7 iDubinski & Carlbera 1991). But, when a dissipa- 
tive gas component is incorporated they tend to be oblate with axis 
ratios of (c/a) Ph = 0.5 and {b/a) Ph > 0.7 jPubinskil 19941) . Fur- 
thermore, haloes formed within a ACDM model are more spher- 
ical with an axis ratio of (c/a) Ph ~ 0.71 while rounder haloes 
with (c/a) 0u ~ 0.77 are obtained for a AWDM model iBullockl 
2002). The shape of our galaxy's halo is disputed and in this re- 
spect the disruption of the Sagittarius dwarf galaxy provides a 
useful test case. By studying tidal tails from this dwarf galaxy 
llbata et all (2001) ruled out the possibility that the halo is signif- 
icantl y oblate (i.e. with an ax is ratio of q Ph < 0.7). However, 
iMartfnez-Delgado et alj J2004I) suggest that the potential of the 
Milky Way dark halo, as constrained by observations of the Sagit- 
tarius dwarf, is likely to be oblate with g$ h = 0.85, while mod- 
elling of radial velocities of M-giants in the Sagittarius dwarf by 
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Figure 3. Luminosity versus mass for dwarf galaxies of the Local Group 
(open circles). The data were taken from Mateo 1 1998) and the luminosi- 
ties were converted to total V-band magnitudes according to Afytot = 
— 2.5 log Ly tot + 4-8. The shaded region shows the mass range of our 
TV-body systems (masses of 2.8 X 10 7 and 5.6 X 10 7 Mq) and the corre- 
sponding possible luminosity range. The luminosity of a satellite of a given 
mass will depend on its mass function and mass to light ratio. 



suggests that the halo is prolate with an average den- 
sity axis ratio within the orbit of Sagittarius close to 5/3. Given all 
this evidence supporting a non-spherical halo, we have adopted an 
oblate one to carry out our numerical simulations with an axis ratio 
of q<s> h = (c/a)* h = 0.8 which corresponds to a flattening in den- 
sity of q Ph ~ 0.73 and q Ph ~ 0.42 in the inner and outer regions, 
respectively. 



3.1.2 Satellite mass models 

Positions and velocities for particles in our self-consistent satellite 
models are randomly drawn from King models (King 1966). The 
latter are a series of truncated isothermal spheres parametrized by 
a concentration defined as c = log (r t /r c ) where r t and r c are the 
tidal and core radii, respectively. These models quite well fit the 
dwarf spheroidal galaxies (Vader & Chaboyer 1994). To completely 
define a satellite we have to specify its mass and tidal radius in 
addition to its concentration (see Table|5j- In all cases, each satellite 
consists of 10 6 equal-mass particles. 

Figure [3] shows the luminosity and mass distribution of real 
Local Group dwarf galaxies. The masses and luminosities (V- 
band) vary between ~ 10 7 and ~ 10 10 M , and ~ 3 x 10 5 and 
~ 10 9 L Q , respectively (see lMatecll99l) . This means that our sim- 
ulated dwarf galaxies fall in the low-mass end of the satellite mass 
distribution. The grey shaded area in Fig.|3|shows the range of pos- 
sible luminosities (given the observational data) for systems with 
masses between 2.8 x 10 7 and 5.6 x 10 7 Mq. The range of pos- 
sible luminosities for our simulated dwarf galaxies is important for 
the discussion in Section l4~3l 
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Table 6. Orbital parameters for our experiments. 



Run 


Satellite 


Pericentre 


Apocentre 


Inclination 




Model 


(kpc) 


(kpc) 


Angle 


1 


Si 


8.75 kpc 


105 kpc 


30° 


2 


Si 


7 kpc 


60 kpc 


45° 


3 
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7 kpc 


80 kpc 


60° 


4 
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40 kpc 


60 kpc 


25° 


5 
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3.5 kpc 


55 kpc 


45° 
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Figure 4. Satellite orbits (see Tablelol projected onto the principal axes of 
our Milky Way galaxy model. The Galactic disc corresponds to the XY 
projection. Distances are in kpc. 



3.2 Numerical tools 

The simulations were evolved with a parallel tree-code for about 
10 Gyr. It uses a tolerance pa rameter of 6 = 0.75 and a fixed inte- 
gration time-step of 1.3 Myr lDubinskill 996). Forces between par- 
ticles are computed including the quadrupole terms and the gravi- 
tational potential uses a Plummer softening length of 35 pc. With 
these parameters the energy conservation in all cases was better 
than 0.1%. All our numerical runs were performed on a Beouwulf 
computer at the IAUNAM-Ensenada. It has 32 Pentium processors 
running at 450 MHz ( Velazquez & Asuilar 2003). Each simulation 
took 96 hours of 'wall-clock' time. 



8 A.G.A. Brown, H.M. Velazquez & L.A. Aguilar 




Figure 5. The distribution of errors in the simulated Gaia astrometric parameters for 1 million stars from the model Galaxy. The grey scale indicates the 
number of stars. Panel (a) shows the parallax errors as a function of V, the spread is caused by the colour dependence of the errors (see Appendix IaI. Panel 
(b) shows the variation of the parallax errors with ecliptic latitude (3 and panel (c) shows the dependence of the errors in on (3. In order to clearly bring 
out the variation over the sky, the magnitude range of the stars was restricted for panels (b) and (c). The discrete steps in (b) and (c) are caused by the discrete 
tabulation of the dependence of the errors on /3. There is no dependence of the errors on ecliptic longitude and <r M „ shows practically no variation with /3. 



3.3 Numerical experiments 

A total of five orbits were chosen to follow up the disruption of the 
satellite models. The parameters that define them are listed in Ta- 
ble |6| The last column refers to the initial inclination angle of the 
orbit with respect to the plane of the Galactic disc. Notice that most 
of these orbits have a pericentre radius inside the Solar radius. Fig- 
ure|4|shows the different projections of these orbits onto the main 
axes of our Milky Way mass model where the plane of the Galac- 
tic disc corresponds to the XY projection. Since the Milky Way 
galaxy is represented by a rigid potential and dynamical friction is 
ignored, the orbital decay of the satellite is not taken into account. 



4 MODELLING THE GAIA SURVEY 

In this section we describe how the stellar positions and velocities 
from the Monte Carlo model of the Galaxy and the TV-body models 
of the satellites are converted to the data that Gaia will deliver: sky- 
positions; parallaxes; proper motions; and radial velocities. We also 
describe how to put together the Galaxy and satellite simulations 
taking care to correctly account for the variation of the number of 
visible stars along the orbit of the satellite, while using as much of 
the iV-body particles as possible. 

4.1 The Galactic data 

To convert phase-space coordinates in the Galaxy model to Gaia 
data we first translate them to the barycentric coordinate frame, 
which has the same axes as the Galactocentric system used for the 
simulations but with its origin at the Sun's position and velocity. 
The barycentric position r b = r— r© and velocity v b = v-V0 are 
converted to the five astrometric parameters and the radial veloc- 
ity Vraxi- The astrometric observables are the Galactic coordinates 
(£, b), the parallax vj, and the proper motions /i£„ = fit cos b and 
Hb- The units of parallax and proper motion are micro-arcsecond 
(/ias) and /ias yr _1 . The radial velocity is given in km s - . 

Subsequently the expected Gaia errors are added to the ob- 
servables according to the accuracy assessment described in Chapt. 
8 of the Gaia Concept and Technology Study Report <ESA| 2000'). 



The dominant variation of the size of the errors is with the apparent 
brightness of the stars in addition to which there is a dependence 
on colour and position on the sky. Thus, as a first step we calculate 
the V-band magnitude of each star from its parallax and assign a 
(V — I) colour based on the spectral type of the star. The relation 
between spectral type and (V — I) was taken from ICo xl 1200(1) . 
Note that no extinction is included in our model Galaxy. 

Like Hipparcos (see lESAlll997t) Gaia will be a continuously 
scanning spacecraft and will accurately measure one-dimensional 
coordinates along great circles in two fields of view simultaneously, 
separated by a well-known angle. The accuracy of the astrometric 
measurements will thus depend on how the parallactic plus proper 
motions of a particular object project onto these great circles. Be- 
cause Gaia will collect its measurements from an orbit in the plane 
of the Ecliptic the astrometric errors will vary over the sky as a 
function of the Ecliptic latitude j3. For the positions and proper 
motions the latter variation of the precision is caused mainly by 
the variation in the number of times an object is observed by Gaia, 
which is dependent on the details of the way Gaia scans the sky 
over the mission lifetime. For the parallax measurements the domi- 
nant factor is the decrease in measurement precision from the eclip- 
tic pole to the equator owing to the way the parallax ellipse of a star 
projects onto a great circle scan (the parallax factor). The precision 
of the radial velocity measurements varies as a function of V and 
spectral type but is assumed not vary over the sky in our simula- 
tions. 

The dependence of the errors on the Ecliptic coordinates 
(A, /3) of a star is taken into account in the simulations by rotat- 
ing the astrometric data to Ecliptic coordinates, adding the errors 
and rotating the result back to the Galactic reference frame. As a 
consequence correlations will be introduced between the errors in 
I and b as well as fit, and /it, even in the absence of any correlation 
in the Ecliptic reference frame. 

Figure [5] illustrates the dependence of the errors in the astro- 
metric data on the apparent brightness of the star and its position 
on the sky. The variation with stellar magnitude is the strongest. We 
show only the variation of parallax error a m with V as the variation 
is similar for the other astrometric parameters. The spread in a m at 
each value of V is caused by the colour dependence of the astro- 
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metric errors (smaller errors for redder stars, see Appendix^). The 
variation of the errors over the sky is shown in Ecliptic coordinates 
(A,/3) and is prominent only for and a" MA „. The errors only 
vary as a function of ecliptic latitude f3. Finally we point out that 
the relative error in the sky positions is at most ~ 10~ 9 (which is 
200 /ias expressed in radians) meaning that distance errors perpen- 
dicular to the line of sight to even the most distant stars (at about 
100 kpc) amount to no more than about 10~ 4 pc. Hence we will 
ignore errors on the sky positions in the rest of the paper. 

Finally, the steps taken to construct the simulated Gaia cata- 
logue are summarized here: 

(1) From r and v calculate the bary centric position and velocity. 

(2) Transform these to the astrometric parameters in the Galactic 
coordinate system and the radial velocity. 

(3) Calculate V from Mv and vj. 

(4) Find the (V — I) colour corresponding to the spectral type 
of the star. 

(5) Rotate the astrometric parameters to the Ecliptic coordinate 
system. 

(6) Add the astrometric errors as a function of V, (V — I) and 
Ecliptic latitude. 

(7) Rotate the resulting astrometric parameters back to the 
Galactic coordinate system. 

(8) Add radial velocity errors. 

For details please see AppendixlAl 
4.2 The satellite data 

Given V and (V — I), the satellite stars are treated in exactly the 
same way as the stars in the Monte Carlo Galaxy where the calcula- 
tion of the astrometric and radial velocity data is concerned. Hence 
here we describe only how the absolute magnitudes and colours of 
the satellite stars are chosen. 

We assume our simulated satellites represent disrupted dwarf 
galaxies and in this and the following sub-section we will make a 
distinction between real dwarf galaxies and our simulated TV-body 
satellites. We will refer to the former as 'dwarf galaxies' or 'satel- 
lites' and to the latter as 'TV-body systems' or 'simulated satellites'. 

All TV-body systems are assumed to represent a single stel- 
lar population of a certain age and metallicity, which forms a 
single isochrone in a colour-magnitude diagram. Combining this 
isochrone with a mass function £(m) for the satellite stars, we as- 
sign to each TV-body particle a value for Mv and (V — I). We use 
the isochrones by Girardi et al. 1 2000). These isochrones have been 
tr ansformed to severa l different photometric systems as described 
in lGirardietalJ<2002h . The procedure for a particular star is to first 
generate its mass and then to interpolate the isochrone in mass to 
find the corresponding values of Mv and (V — I). 

The colour-magnitude diagrams for a number of Local Group 
dwarf galaxie s are very similar to those of Galactic globular clus- 
ters ( see e.g., iFeltzing. Gilmore & Wisdll999t iMonkiewicz et alJ 
1999). This suggests that we can use a globular cluster like mass 
function for our TV-body systems. Based on the study of a dozen 
globular clusters lParesce &De Marchill2000h propose a lognormal 
mass function for these systems for stars below about 1 Mq . This 
mass function has a characteristic mass of 0.33 Mq and a standard 
deviation of 0.34. The slope of the mass function between ~ 0.3 
and ~ 1 Mq is similar to that of a power-law mass function for 
which f (m) oc m~ a , with a ~ 1.5. For N-body snapshots older 
than 5 Gyr the stars in the simulated satellites will all have masses 
less than about 2 Mq and for simplicity we use a single power-law 



mass function £(m) oc m -1,5 . Hence we may be overestimating 
the number of low-mass stars (below 0.3 Mq) in our TV-body sys- 
tems with respect to the more massive stars. 



4.3 Converting satellite TV-body particles to stars 

Having made a considerable effort to realistically — in terms of the 
number of stars — simulate the Galaxy as it will be seen by Gaia, 
we want to ensure that the satellite TV-body particles are properly 
converted to simulated stars representing a dwarf galaxy. Getting 
this right is not trivial and we explain below why this is the case 
and describe in detail our solution to this problem. 

Given a certain distribution of stars along the orbit of a par- 
ticular dwarf galaxy (corresponding to an TV-body 'snap-shot'), the 
number of stars from this satellite that will end up in the Gaia cat- 
alogue depends on three factors: 

(i) The overall number of stars (i.e., luminous particles) in the 
dwarf galaxy. This number is determined by its overall luminosity 
and the stellar mass function. 

(ii) The Gaia survey limit. This magnitude limit determines the 
faintest dwarf galaxy stars visible by Gaia. 

(iii) The variation of the number of visible satellite stars along 
the orbit. This variation is caused by the variation in distance from 
the dwarf galaxy stars to the Sun. 

If we take the TV-body particles as representing a sampling of 
the distribution of the real satellite stars along the orbit, then the 
expected number of satellite stars visible in the Gaia survey can 
be calculated as follows. The overall number of stars (luminous 
particles) TVj un , in a real dwarf galaxy can be obtained from its 
total luminosity in the V-band Lv.tot and its mass function £(m): 



TVlun 



(Lv) 



(Lv) 



J Lv(m)(,(m)dm 



(8) 



Lv (m) is the V-band luminosity of a dwarf galaxy star with mass 
m and (Lv) is the average stellar luminosity in the dwarf galaxy. 
At each distance s; at which an TV-body particle is located, the frac- 
tion of visible satellite stars is set by the magnitude limit Mv,max 
which is given by: 



Mv,max(Si) = Vum - 5 log S t + 5 , 



(9) 



where Vn m is the Gaia survey limit. The fraction of visible satellite 
stars fi at distance Si is then: 



fi(si) = 



/m u - 
m(i 



™(AV,max) 



£(m)dm 



(10) 



where m up is the upper mass limit of the mass function. The overall 
fraction / of visible satellite stars is then given by: 



/ = 



TV sim 



(11) 



where TV s i m is now the total number of simulated satellite particles. 
Hence the number of satellite stars expected in the Gaia survey 
is /TVi um- Note that we assume that the mass function does not 
vary along the orbit of the satellite. This may not be realistic if the 
satellite preferentially looses lower mass stars owing to tidal forces. 

Now, for a normal mass function the number of low-mass stars 
will heavily dominate, which means that the fractions fi can be 
very small depending on the distance Si. Figure |6| shows an ex- 
ample of this for a snapshot at about 10 Gyr of TV-body run no. 
1 (see Table|6j. The stellar population of the simulated satellite is 
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distance (kpc) 

Figure 6. Visible fraction of the simulated satellite stars at each TV-body 
particle position as a function of the distance to the observer for TV-body 
run no. 1. The shaded histogram (vertical scale on the right) shows the dis- 
tribution of distances Sj. The solid line (vertical scale on the left) shows 
for each distance Si the fraction of visible stars /; calculated as in Sec- 
tion !4.3l The features in this line reflect features in the luminosity function. 
The dashed line shows the fraction of visible stars when taking into account 
only the stars brighter than the faintest satellite star that enters the Gaia sur- 
vey. The fractions of visible stars were calculated using a 10 Gyr isochrone 
of metallicity Z = 0.004 and a mass function oc m -1,5 . The overall 

fraction of visible TV-body particles is only 2.5 X 10~ 3 and 10~ 2 for the 
two cases considered. 



represented by a 10 Gyr isochrone for a metallicity Z = 0.004, 
with a mass function £(m) oc m~ 1,5 . Three things are shown as a 
function of the distance Si of each TV-body particle i. The shaded 
histogram shows the distribution of distances S; and illustrates the 
wide variety of distances at which the simulated satellite particles 
are located. The solid line shows for each distance s, the fraction of 
visible stars fi calculated as outlined above. The TV-body particles 
are taken to represent a fully populated mass function which in this 
case ranges between 0.15 and 0.88 Mq. Note that fi is nowhere 
larger than ~ 0.1, while the overall fraction / of visible stars is 
2.5 x 10 -3 . This means that if we assume that our TV-body parti- 
cles represent stars over the full range of the satellite mass function, 
only 2 500 out of the 1 million simulated stars will end up in our 
synthetic Gaia catalogue! 

This is clearly a very wasteful way to make use of the (compu- 
tationally) expensive TV-body calculations and the question is how 
to retain more TV-body particles in the synthetic Gaia catalogue 
whilst at the same time preserving the variation of the visible frac- 
tion of stars along the TV-body system's orbit. One straightforward 
improvement can be made by realizing that for mass functions that 
extend all the way to late M-dwarfs most of the generated synthetic 
satellite stars are too faint to ever make it into the Gaia survey, re- 
gardless of the position at which they are located along the satellite 
orbit. To eliminate this problem we can first find the position in the 
TV-body snapshot that is closest to the Sun (at a distance s m i n ) and 
calculate the corresponding absolute magnitude, My.faintcst, of the 
faintest satellite stars visible by Gaia. Subsequently in assigning 
luminosities to the simulated satellite stars we use the luminosity 
function only for My ^ My, f a i nt0 st • The result is shown in Fig.|S| 




55 60 
distance (kpc) 

Figure 7. Same as Fig.lolfor the simulated satellite from TV-body run no. 4. 
A 5 Gyr isochrone and the same mass function were used. Note the much 
smaller range in distances for this TV-body system, resulting in a much 
larger fraction of visible stars if the luminosity function is cut off at My ^ 

M v 

faintest - The overall fraction of visible stars is 3.4 x 10 ^ when using 
the full mass function and 0.78 when only stars with My ^ TV/\/ faintest 
are considered. 



by the dashed line. The visible fraction of stars increases every- 
where by more than an order of magnitude due to the exclusion 
of the faintest (never visible) stars. However the overall fraction of 
visible stars is still very low at 9.8 x 10~ 3 . This is a consequence of 
the large variation in distances for the orbit of the simulated satellite 
from TV-body run no. 1. Figure^shows the distribution of distances 
and values of fi for the satellite from run no. 4 at 5 Gyr (using a 
5 Gyr isochrone and the same mass function), which has a much 
more circular orbit with much less distance variation (compare the 
orbits in Fig.|4j- For this TV-body system one can retain about 78 
per cent of the TV-body particles if the luminosity function cut-off 
at M v ^ My.faintcst is used. 

A further improvement is now obvious. Instead of cutting off 
the luminosity function at My.famtcst we choose a brighter limit 
and retain more of the TV-body particles in the synthetic Gaia cata- 
logue. This can be interpreted as simulating the fact that in practice 
one will use bright 'tracer' stars (such as those at the tip of the 
AGB or RGB) to map out the satellite in phase-space. That is, we 
assume that the TV-body particles only represent stars at the high 
mass end of the satellite's mass function. In fact most of the present 
surveys for substructure in the Galactic halo rely on using tracer 
populations that are more easily isolated photometrically in colour- 
magnitude diagrams (see e.g. , iMajewsk^t alll2003l iMartin et alJ 
l2004tlRocha-Pinto et all2004tlHelmj|2004l) . 

In principle the luminosity function cut-off magnitude 
My^raccr can be chosen bright enough to retain all TV-body par- 
ticles. However, there are two constraints on its value, the first of 
which is the obvious limit set by the magnitude of the brightest star 
in the luminosity function. The second constraint is set by the desire 
for a realistic simulation of a satellite as it will appear in the Gaia 
catalogue. By assuming that the TV-body particles represent only 
a tracer population we are effectively simulating a more massive 
dwarf galaxy when converting the TV-body results to Gaia obser- 
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Figure 8. Visible number of stars as a function of the luminosity function 
cut-off My tracer f° r a satellite with its stars distributed along the orbit 
of TV-body run no. 1 . The dashed curves with triangles show the numbers 
for 'real' dwarf galaxies (i.e., with fully populated luminosity functions) 
of total luminosities: My tot = — 17, —15, —13, —11 (where My tot = 
— 2.5 log Ly.tot + 4.8). The dotted line shows the value of My faintest 
beyond which the dashed curves become flat. The solid curve with dots 
shows the numbers of visible stars for the ./V-body system, where now all 
stars are assumed to be brighter than My tracer- The fractions of visible 
stars were calculated using a 10 Gyr isochrone of metallicity Z = 0.004 
and a mass function £( m ) <x m . 



vations. In reality we have simulated rather small dwarf galaxies 
(see Fig. |3 which cannot be scaled up trivially to more massive 
systems because the latter respond differently to the Galactic tidal 
forces. Thus a balancing between the number of retained TV-body 
particles and the total mass (or luminosity) of a real dwarf galaxy 
that we wish to simulate is required. 

This balancing is illustrated in Figs. [8] and |5| for the TV-body 
systems discussed above. The figures show for the TV-body system 
and for the corresponding real dwarf galaxies how the number of 
stars visible by Gaia changes as the magnitude cut-off of the lumi- 
nosity function, My tracer, changes. These numbers are calculated 
according to equations l8t — 1 1 1 i - For a real dwarf galaxy the number 
of visible stars will increase as the value of Mytracer goes up, re- 
flecting a probing of the luminosity function to fainter magnitudes. 
At some point My,tracer becomes larger than My,fajntest and the 
number of visible stars stays constant (probing deeper into the lu- 
minosity function only adds stars invisible to Gaia). This is shown 
by the dashed curves and triangles, where each curve represents a 
dwarf galaxy of different Lytot (or total mass). The dotted line 
shows the value of My,f a mtest beyond which the curves are flat. 

For the TV-body system the number of visible particles will 
decrease as the value of My,tracei goes up. This is because now all 
particles are assumed to be brighter than My,tracer which leads to 
a larger fraction of simulated stars that are fainter than the survey 
limit as the luminosity function cut-off becomes fainter. In Figs. [8] 
and [9| the solid curve with dots represents the visible number of 
stars for the TV-body system. Note that this curve continues to drop 




Figure 9. Same as Fig.|5]for a satellite with its stars distributed along the 
orbit of from ./V-body run no. 4. A 5 Gyr isochrone and the same mass 
function were used. Note the brighter value of My faintest which reflects 
the larger average distance of this satellite. 



beyond My,tracer = Myfamtcst which reflects the increasing num- 
ber of simulated stars that will never enter the Gaia survey. 

The balancing of the number of retained TV-body particles and 
the total mass of the dwarf galaxy we wish to simulate thus con- 
sists of choosing a dwarf galaxy luminosity and then finding the 
value of My,traccr for which the number of expected tracer stars 
from this satellite in the Gaia survey is equal to the number of re- 
tained TV-body particles. For TV-body run no. 1 one can simulate 
an Mytot = —17 satellite by choosing Mytracer ~ 0.5, which 
then results in between 7 x 10° and 8 x 10 TV-body particles en- 
tering the synthetic Gaia catalogue. Note that in practice one can 
choose a value of My,tracer corresponding to an approximate value 
of Afytot by inspecting diagrams such as Figs. [8]and|9| and then 
calculate the actual value of My jto t based on the amount of TV- 
body (tracer) particles that enter the simulated Gaia survey. 

We now return to Fig. [3] which shows that for the mass range 
of the satellites we have simulated the total V-band magnitude can 
be at most My, to t ~ —15. From Figs. [8] and [5] one can see that 
this luminosity is simulated for My,tracet ~ 3 and that then we 
retain about 2 x 10 5 TV-body particles for each satellite. We still 
consider this wasteful and decided to use My,tracer = 0.5 in order 
to retain many more TV-body particles. This implies an unrealis- 
tic value of My ~ —17 for the absolute magnitude of the sim- 
ulated satellites in the V-band. The main problem being that the 
corresponding masses of more than 10 9 Mq (see Fig. |3J would 
require simulations with at least an order of magnitude more parti- 
cles and with the effects of dynamical friction included. This was 
beyond our computational capabilities at the time the simulations 
were done. However, the purpose of this paper is not to provide 
detailed and completely realistic simulations of the Galaxy and its 
satellites. The main purpose is to get an insight into the challenges 
one will have to deal with when searching for the debris of past 
merger events in the Gaia data. Therefore we feel it is more impor- 
tant to get a good sampling of the distribution of satellite stars in 
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Figure 10. Mollweide projections of the distribution on the sky of the stars 
in our simulated Gaia catalogue. The top image shows only the satellites 
projected on the sky (from simulation runs no. 1 and 5 at 10 Gyr, where an 
extra satellite was generated by reversing the angular momentum for run 1). 
The bottom image shows the Galaxy and satellite stars combined. Note how 
the satellite stars become almost totally invisible. The only conspicuous 
over-density is in the top left corner of the image. The dark rectangle is the 
area on the sky with coordinates -90° < I < +90° and -5° < b ^ +5° 
that was left out of the simulations . The figures were generated using the 
Healpix iGorski. Hivon & Wandelt 1999) division of the celestial sphere 
into cells of equal area and counting the stars inside each cell. Dark areas 
correspond to empty cells while the individual dots indicate cells with 1 
or more stars. Areas where all contiguous cells are filled appear as a grey 
scale. 

the Gaia catalogue, by artificially increasing their numbers, than to 
construct a completely realistic simulation of the satellite debris. 

Finally, we summarise here the steps for simulating the satel- 
lite galaxy as it will appear in the Gaia catalogue: 

(1) calculate the distance to each TV-body particle. 

(2) Choose a value My.tracer for the minimum brightness of the 
tracer stars represented by the TV-body particles. 

(3) Re-normalise the luminosity function between My(m up ) 
and A/v.traccr (i.e., all TV-body particles represent stars brighter 

than Mv,tracer). 

(4) Generate simulated stars for each TV-body particle position 
(distance) according to the renormalized luminosity function and 
decided whether or not they are bright enough to enter the Gaia 
survey. 

(5) For stars bright enough to enter the Gaia survey assign a 
(V — I) colour by interpolating in mass along the isochrone. 

(6) Generate the Gaia data according to the steps outlined in 
Sectionl4~T1 

This procedure ensures a properly sampled luminosity function 
along the orbit of the dwarf galaxy. 

Our procedure assumes that one can select the bright tracer 



stars a priori which may not always be possible in practice and may 
not be desirable. A way to get around this problem and still use as 
many TV-body particles as possible is the following. Assume that 
all TV-body particles are brighter than T\i"v,famtcst and calculate the 
predicted number of TV-body particles that make it into the survey 
as for Figs. [8j and |9| Calculate the ratio R = TV a i m /TV v is, where 
TVvis is the number of visible particles for Mv,traccr = TWv,f a i n test ■ 
Subsequently repeat steps (l)-(6) above R times for each TV-body 
particle. That is, generate R stars for each TV-body particle and re- 
tain those that are bright enough to enter the survey. This way all 
10 6 TV-body particles can be retained in the simulated Gaia data 
whilst still retaining a proper luminosity function sampling along 
the orbit of the dwarf galaxy. However one then faces the prob- 
lem of possibly having simulated a system that is too massive (but 
this can be tuned by not retaining all 10 6 particles) and of having 
simulated R stars at exactly the same position with the same ve- 
locity vector. It is not clear how one would redistribute these extra 
stars along the dwarf galaxy orbit such that their energy and angu- 
lar momentum are consistent with that of the TV-body particles. We 
therefore decided not to pursue this option further. 

To close this section we show in Fig. llOl the sky distribution 
of the stars in our simulated Gaia catalogue. The simulated satel- 
lites included in these figures correspond to TV-body runs 1 and 5 
of Table|6|where a third satellite was generated by reversing the an- 
gular momentum of the satellite from run 1 . The main thing to note 
is how the satellite stars are completely swamped by the Galactic 
stars in these sky projections. This is not a new fact but serves to 
illustrate the need for clever search techniques to find the debris of 
disrupted satellites, an issue to which we turn in the next section. 



5 RESULTS FROM GAIA SURVEY MODELLING 

iHelmi & White! 1 1999ft showed that it will be quite difficult to find 
spatial correlations in the debris originating from satellites, whose 
disruption occurred in the early evolution of the Milky Way after 
several passages through pericentre. In this respect, techniques to 
identi fy streamers lik e the 'Great Circle Cell Counts' proposed by 
Ijohnston et alj 1 19961) will be of little practical use. This method 
requires that the orbital plane of the satellite remains almost un- 
changed during its disruption process, so that the debris occupies 
great circles on the sky as seen from the Galactic Centre. This only 
occurs if the outer halo is spherical or if the satellite follows a polar 
or planar orbit in t he cas e of an oblate halo. 

IHelmi et alJ 1 19991) showed that useful information can be 
obtained from the phase-space structure of these fossil remnants 
which allows tracing the merger history of our Galaxy. Initially 
satellites are clumped in both the configuration and velocity space, 
and also in the integrals of motion space. The method of search- 
ing for debris in the integrals of motion space relies on the ba- 
sic assumption that these quantities are conserved, or evolve only 
slightly, and hence their lumpiness should be preserved. For a 
spherically symmetric halo the integrals of motions are the energy, 
total angular momentum and angular momentum around the Galac- 
tic Z-axis, (E, L, L z ); for an axisymmetric halo E and L z are 
preserved while L evolves owing to the orbital precessing of the 
satellite but it still exhibits some deg ree of coherence. The integr als 
of motion technique was applied bv lHelmi & de Zeeuwl ^2000) to 
an artificial Gaia halo catalogue generated from TV-body simula- 
tions of the disruption of satellites in a Galactic potential. From this 
study, they conclude that identifying substructure within the stellar 
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Figure 11. Error propagation in the E-L z plane. The effect of the astromet- 
ric and radial velocity errors on the observed value of E and L z is shown 
using the example of a particular star with r and v corresponding to the 
satellite from simulation run 1 at 10 Gyr. Each line indicates how E and 
L z change from the true value (indicated by the dot) as the observables 
are varied within ±3<r from their true values. The values of the astromet- 
ric parameters, the radial velocity and the corresponding errors are listed in 
the legend. The solid line shows the effect of the error in the parallax. The 
dotted line shows the effect of errors in fi£ t . The dashed line indicates the 
effect of radial velocity errors. The inset enlarges the region around the true 
value of (/'<', Lz) in order to better show the effects of the proper motion 
and radial velocity errors. Note that the error in /i b has almost no effect 
on (E, L z ) for this particular star. The scale of the figure is the same as 
for Fig. 1121 For convenience the potential energy was computed using the 
potential from equations (8)-(10) in Helm i & de Zeeuwl l2OO0t) . 



halo resulting from merger events in the past will be relatively easy 
even after Gaia observational uncertainties are taken into account. 

However lHeimi & de Zeeuwl <2000 ) made two important sim- 
plifying assumptions. They constructed a synthetic Gaia catalogue 
consisting of stars that are intrinsically as well as apparently bright 
(My = 1, V ^ 15), and which are all located relatively near 
the Sun (d < 6 kpc). The limit on apparent brightness was set 
by the demand of having good radial velocities available for their 
sample. Secondly, their simulated halo consists exclusively of in- 
falling satellite Galaxies. These two assumptions lead to a sam- 
ple of halo stars which have small observational errors and will be 
highly clumped in the integrals of motion space, making the disen- 
tangling of the different debris streams relatively easy. 

In reality there may well be a 'background' population of halo 
stars that have a smooth distribution in (E,L,L Z ) space owing 
to previous phase mixing and dynamical friction or because they 
constitute a component that was formed in a monolithic collapse 
and thoroughly mixed by violent relaxation. In addition real de- 
bris streams will have a spread in stellar luminosities (even if only 
bright tracers are selected such as in our simulations) and may be 
located much further away. These will lead to larger observational 
errors which will propagate into (E, L, L z ) space. Both effects will 
smear out the signatures of individual debris streams. We demon- 



5.1 Error propagation in E-L z space 

Figure ITU schematically illustrates how the Gaia observational er- 
rors propagate into the E-L z plane. As an example a star from the 
satellite of simulation run no. 1 is chosen and the parallax, proper 
motions and radial velocities are varied by ±3<j from their true val- 
ues. The lines in the diagram trace the effect of varying each of the 
observables individually. The measured parallax can in principle 
become negative but then the calculation of the observed angular 
momentum becomes ill-defined. Hence Fig. II ll onlv shows the ef- 
fect of the positive branch of the measured parallaxes. The paral- 
lax error is by far the dominant contribution to the errors in E-L z 
space and the observed values of (E, L z ) are smeared in a preferred 
direction. In addition it can be seen that the angular momentum can 
change its sign for large parallax errors. This happens because the 
star, which is located beyond the Solar circle in the direction of 
I — 320°, can get placed inside the Solar circle if the measured 
parallax is much larger than the real parallax (measured distance 
too small). This puts the star at an observed position on the oppo- 
site side of the Galactic centre which for an unchanged direction of 
its observed velocity leads to a reversed angular momentum. 

Another effect is caused by the non-linear dependence of E 
and L z on parallax, i.e. both quantities depend on functions of 
1/zn (distance). This leads to systematic errors in the sense that 
the expectation value for the observed E or L z does not equal 
the true values of these quantities. This is a well known problem 
when dealing with the determination of dista nces to stars using 
parallaxes with relat ively large errors (see e.g.. lBrown et alJl997t 
ISchr6deretaiT2 004). The detailed systematic effects in the E-L z 
plane depend on the actual distribution of the observed stars in 
space, on the distribution of the parallax errors and, importantly, 
on the way the sample of stars is selected. Any selection made on 
the size of the parallax or on its quality (a-^jvj) will skew the un- 
derlying distribution of true distances which will be reflected in 
the calculated integrals of motion. A detailed investigation of these 
systematics is beyond the scope of this paper but is important in 
the future use of Gaia data to study the integrals of motion in the 
Galaxy. Appendix|B|contains the equations for calculating the an- 
gular momentum of a star from astrometric and radial velocity data. 
They can be used to further study the effects of measurement errors. 



5.2 Searching for satellite remnants in E-L z space 

We constructed E-L z diagrams for all of our satellite TV-body sim- 
ulations together with the synthetic catalogue of the background 
population associated with the Milky Way. For clarity in the follow- 
ing discussion we concentrate on just a few satellites. The diagrams 
are shown in Fig. 1121 The upper panels refer to the E-L z diagrams 
without Gaia errors (from left to right: Galaxy only, satellites only 
and Galaxy and satellites combined). The middle panels show the 
effect of adding the expected Gaia observational errors. In this case 
the sample was restricted to stars with positive parallaxes and mea- 
sured radial velocities in order to calculate E and L z . This selec- 
tion criterion introduces an implicit limit on the apparent bright- 
ness of the stars of V < 17-18. In the bottom panels the sample of 
stars selected from the catalogue is restricted to stars brighter than 
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Figure 12. E-L z space for three of our disrupted TV-body satellites and the synthetic Milky Way catalogue. The upper panels show the error free E—L z 
diagrams (from left to right: Galaxy stars only, satellite stars only, and Galaxy and satellite stars combined, respectively). The middle panels show E-L z space 
after the expected Gaia observational errors have been added. Only stars with radial velocities and positive measured parallaxes are plotted. Finally the lower 
panels show only the simulated stars from a high-quality sample for which V < 15 and to / a m > 5. These figures have been constructed as two-dimensional 
histograms by counting stars in 250x250 cells that divide E-L z space. The white areas correspond to empty cells while the individual dots indicate cells with 
1 or more stars. Areas where all contiguous cells are filled appear as a grey scale. 



V = 15 having a parallax signal-to-noise of w ja^, > 5 . In these 
diagrams, the entire Galaxy catalogue consists of 3.5 x 10 8 stars 
and the satellite models correspond to runs 1 and 5 of Table|6|where 
a third satellite was generated by reversing the angular momentum 
of the satellite from run 1. All satellites have Afv.traccr = 0.5 (i.e. 
Afv.tot ~ —17) and the stellar luminosities and colours were com- 
puted using a 10 Gyr isochrone with metallicity Z = 0.004 and a 



mass function ^(m) oc m . Energies were computed employing 
the Galaxy model described in section|3| 

In practice one will not know the actual potential of the Milky 
Way and thus an estimate of this quantity is required to calculate 
E. To simulate the effect of an erroneous estimate of the potential 
we also computed the energy using an alternative Galaxy potential 
model (see equations (8)— (10) of Helmi & de Zeeuw 2000). This 
results in a larger spread of the Galaxy and satellite stars in energy 
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Figure 13. Energy histograms obtained from E-L z space for all satellite and Galaxy stars by taking slices in the integral of motion L z . The leftmost two 
/^-histograms correspond to all stars located in the interval —7 X 10 3 < L z < containing the prograde satellites (from left to right the entire and restricted 
Gaia samples, respectively). Similarly, the right two B-histograms were constructed by using the interval 10 3 < L z < 7 X 10 3 enclosing the retrograde 
satellite. In all cases, thin lines indicate error-free data while thick lines correspond to data after Gaia errors are added. Prograde satellites are centred at 
energies E fs —6.5 X 10 4 (km s -1 ) 2 , E f» —4.2 X 10 4 (km s -1 ) 2 ; the retrograde satellite is centred at energy E —4.2 X 10 4 (km s -1 ) 2 . 



and this spread can be used in combination with the angular mo- 
mentum measurements and a detailed modelling of the identified 
debris streams in order to constrain the Galactic potential. How- 
ever the E-L z diagrams with the observational errors included 
look largely the same as those in Fig. ^| so we will not consider 
this issue further in the discussion that follows. 

It can be seen that satellite remnants are easy to identify if 
there are no observational errors, even in the presence of a stellar 
background (upper panels of Fig. I12i . However, once the obser- 
vational errors are introduced, the identification process becomes 
more difficult. Notice how the satellite and Galaxy stars are now 
smeared out in preferential directions in the E-L z diagram, which 
is caused by the dominance of the parallax errors as illustrated in 
Fig. 1111 The middle panels of Fig. ^| demonstrate that if the en- 
tire Gaia Survey is taken into account it becomes very difficult to 
disentangle the different merger events. In particular, satellites on 
prograde orbits become buried in the very crowded region occupied 
by the Galactic disc. The identification process can be significantly 
improved when the search in E-L z is restricted to a 'high quality' 
sample as illustrated in the lower panels of Fig. 1121 In this case, 
even a visual identification clearly shows all three satellites despite 
the presence of the Galactic stellar background. In this respect, the 
middle dia gram of the lower panels in Fig. 1 1 2l is similar to Fig. 4 
of lHelmi & de Zeeuwl 1 2000 ). However the clumping of the debris 
streams in E-L z is less clear and they are connected to each other 
by the background population of Milky Way stars. This makes the 
application of a 'friends of friends' algorithm to search for debris 
streams (as proposed by Helmi & de Zeeuwl2000l) less attractive. 

Another very important feature of these E-L z diagrams is that 
the observational errors introduce a lot of false high density struc- 
tures associated with the satellites (see middle column in Fig. ll2i . 
These structures look like singularities or caustics in phase-space 
and from a naive analysis could be mistaken for features in the 
topology of phase-space corresponding to physical entities or past 
merger events that in reality do not exist. 

Taking slices in the integrals of motion by selecting certain 
ranges of E or L z allows us to make a more quantitative assess- 
ment by focusing our attention on E or L z -histograms containing 
the relics of the disrupted satellites. A similar approach has been 
used by M eza et all 12005) for the case of the debris of the wCen 
dwarf, the galaxy that presumably contributed the uo Cen globular 
cluster to our galaxy. They were able to identify histogram peaks 
coinciding with the loss of satellite particles in the last three peri- 
centric passages (see their fig. 4). However, our case proves to be 
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Figure 14. L z histograms obtained from E-L z space for all satellite and 
Galaxy stars by taking a slice in the integral of motion E. The left panel 
shows the L z -histogram considering the whole survey in the energy interval 
— 10 s < E < while the right panel shows the distribution of L z for 
our high quality sample in the same energy range. The lines have the same 
meaning as in Fig. 1131 



more challenging, as can be appreciated in Figs. ll3l and ll4l Here 
the thin lines refer to error-free data and the thick lines to data after 
Gaia errors have been taken into account. E histograms for error- 
free data show a double peak around each satellite. These peaks 
seems to be related to stars coming from the leading and trailing 
tails of the disrupted satellite. However, these peaks are smeared 
out over the entire energy space when Gaia errors are added, while 
only some bumps survive in our high quality sample (notice that 
even in this more favourable situation any trace of one of the pro- 
grade satellites has been completely erased). A similar behaviour is 
seen in the L z histograms, although in this case all three satellites 
are preserved after the errors are added in the high quality sample. 

Finally, we note that dynamical friction has been ignored in 
our numerical simulations of satellite disruption. If dynamical fric- 
tion is taken into account L z will not be conserved and its evolu- 
tion will strongly depend on the initial satellite mass as well as on 
the orbital angular momentum with respect to the angular momen- 
tum of the main Gal axy. For exam ple prograde orbits decay faster 
than r etrograde ones (Huang & Carlberg 1997; Velazquez & White 
1999). However, Helmi, White & Soringel (2003) have pointed out 
that the conservation of phase-space density remains a promising 
technique for providing useful information about the Galaxy for- 
mation history even in a more realistic situation where an entire 
galaxy builds up from mergers in a A-CDM cosmology. 
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Figure 15. Radial velocity versus distance diagrams for our simulated Gaia catalogue. The same satellites as in Fig. I12l were used and superposed on the 
Galactic background. The left panels shows the data without observational errors. The middle panel shows the effect of adding the Gaia measurement errors, 
and the right panel shows the Gaia data for a high quality subsample restricted in magnitude, V < 15, and parallax, signal-to-noise ro/<r ro > 5 as in Fig. 1121 
Note the change in the horizontal scale in this last panel. These figures were generated in the same way as Fig. 1121 



5.3 Searching for satellite remnants in r-v ra d space 

Johnston. SDersel & Hernauist 1 1995) suggested that diagrams of 
radial velocity versus distance would provide a powerful tool to 
look for fossil signatures in the halo of our Galaxy. These r-u ra d 
diagrams are related to the E-L z diagrams since, roughly speak- 
ing, the total energy of particles reflects the apocentre of the orbit 
while the pericentric radius is associated with the orbital angular 
momentum. This technique has not been applied in practice ow- 
ing to a lack of data sets with enough radial velocities and reliable 
distances. 

Figure [T5l shows the r-v r ^ diagrams for our simulated Gaia 
catalogue. The same satellites were used as for the E-L z diagrams. 
The data without the observational errors are shown in the left 
panel. The debris of the disrupted satellites can clearly be distin- 
guished as narrow structures in this diagram but it is not possible to 
associate the latter with one or another of the satellites. When the 
observational errors are included most of these fossil signatures are 
washed out and a straightforward detection will be quite difficult 
as can be appreciated in the middle r-f ra d diagram. Restricting the 
data to the same high quality subsample as defined above, does not 
change the situation although it is possible to discern some struc- 
tures related to the tidal debris left by the disrupted satellites. 

5.4 The need for complementary search methods 

These results show that a straightforward search for the debris of 
satellite galaxies in E-L z or r-u ra d space may not be very effi- 
cient in practice. A good start can be made by limiting the search 
to high quality samples, where the contrast between satellites and 
Galaxy can be further enhanced by concentrating on certain regions 
of the sky (away from the Galactic disc for example). However the 
effects of error propagation and selection biases (see Section fSTTt 
should be carefully accounted for. Once a first set of stars belong- 
ing to a particular remnant has been identified in E-L z one can in 
principle trace the rest of the stream by modelling the orbit of the 
satellite and by making use of the astrophysical properties of the 
stars (age, metallicity, alpha-element enhancements) which will be 
known from the photometric data provided by Gaia. Alternatively 
one can select samples of metal poor stars (thereby also lowering 
the Galaxy-satellite contrast) from the Gaia catalogue in order to 



study the very oldest remnants. Gaia will also provide extremely 
well calibrated photometric distances based on the parallaxes of 
the nearer and/or brighter stars. These calibrations can be used to 
obtain accurate distances for faint stars for which Gaia parallaxes 
are of poor quality. Better distances for the fainter stars will dramat- 
ically improve the quality of the E~L Z and r-i> ra d diagrams. These 
complementary procedures can be studied with our simulations by 
adding a more complete model of the Gaia photometric data. 



6 CONCLUSIONS 

Gaia will provide us with an unprecedented wealth of information 
on the structure and formation history of the Milky Way galaxy, for 
the first time providing full phase-space information across its en- 
tire volume. Exploiting this data set for studies of the assembly of 
the halo and other components of the Milky Way will be challeng- 
ing owing to the enormous volume of data; astrometry and pho- 
tometry for 1 billion stars and radial velocities for 100-150 million 
stars, over the entire sky. A good preparation will require familiar- 
izing ourselves with the handling of such large datasets as well as 
exploring the best ways of disentangling remnants of past merger 
events from each other and from the Galactic background popula- 
tion. 

We have started to address this problem by simulating the 
Gaia catalogue with a realistic number (10 8 -10 9 ) of entries. This 
was done by combining a Monte Carlo model of the Milky Way 
with A-body simulations of satellites being disrupted as they or- 
bit in the potential of our Galaxy. The Monte Carlo model of the 
Milky Way consists of 3.5 x 10 8 stars and was generated by taking 
into account how the Gaia survey limit introduces magnitude lim- 
ited spheres for each stellar type. The region on the sky limited by 
Galactic coordinates -90° +90° and -5° b ^ +5° was 

left out of the simulations for practical (computation time) reasons. 

The models of the debris streams were generated by simulat- 
ing the disruption of dwarf galaxies orbiting our galaxy. The Milky 
Way is represented by a rigid potential and a tree-code was used 
to carry out the A-body simulations of the orbiting satellites over 
a time period of about 10 Gyr. The simulations were carried out 
for satellites placed on five different orbits which vary in apocen- 
tre, pericentre and the initial inclination with respect to the Milky 
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Way's disc. Multiple debris streams can then be simulated by com- 
bining different TV-body snapshots (different orbits and/or ages), 
and each snapshot can be rotated around the Z-axis, flipped with 
respect to the disc plane, or reversed in angular momentum. 

We describe in detail how to correctly combine the Milky Way 
and satellite models such that the variation of the number of visible 
stars along the orbit of a satellite is taken into account. This is most 
easily achieved by assuming that the TV-body particles represent a 
bright tracer population of a larger underlying dwarf galaxy. 

After combining the models the synthetic Gaia catalogue was 
generated by adding the astrometric and radial velocity errors ac- 
cording to the prescriptions resulting from detailed preparatory 
studies of the Gaia mission. The observational errors were mod- 
elled taking into account the dependence on apparent magnitude, 
colour and sky position of the stars. To assign photometric prop- 
erties to the particles we use a Hess diagram for the Solar neigh- 
bourhood for Galactic particles, while for the satellite particles we 
use isochrones from the Padova group. Our work thus introduces a 
thorough methodology to incorporate TV-body simulations within a 
simulation of an actual observational survey in a manner that ad- 
dresses questions of kinematics, metallicity and photometric prop- 
erties of the simulated system, as well as errors in the observables. 

The catalogue was used to study the appearance of the Galaxy 
and satellite data in E-L z space. The results show that it will be 
challenging to trace the debris streams throughout the entire vol- 
ume of the Galaxy when Gaia errors are properly accounted for. 
However, limiting the sample to bright stars with good parallaxes 
allows a preliminary search to be done which should be followed 
up by complementing searches in E-L z by, for example, a further 
tracing of the debris streams based on the astrophysical properties 
of the stars which will be known from the photometry provided by 
Gaia. In addition, using photometric distances, which will be very 
well calibrated after the Gaia mission, will yield accurate distances 
for faint stars which can be used instead of the parallaxes, thus 
avoiding the problems associated with the propagation of parallax 
errors into the integrals of motion. 

Ultimately our aim is to use these simulations to study much 
more thoroughly and realistically the possible search techniques 
one can use to recover the remnants of satellites in the stereoscopic 
and multi-colour data set that Gaia will provide. This includes a 
more complete study of the use of the integrals of motion method 
by also using the full angular momentum vector L and a further 
exploration of the radial velocity versus distance diagrams. In ad- 
dition we will explore how the astrophysical information from the 
Gaia photometry can be used in the search for debris streams. This 
requires a number of improvements to our simulations, which in- 
clude: 

• A proper model of the Gaia photometric system which is 
available using the tools developed by the Gaia Photometric Work- 
ing Group (see Jordi & H0 j2005tlCarrasco et all2005l) . 

• Updating the predicted Gaia observational errors to the lat- 
est available assessments. Specifically, the radial velocity errors we 
used are some what optimistic and a better assessment of them is 
now available jKatz et alfeool . 

• Using a clumped distribution of stars in the halo. The simu- 
lation of a clumpy halo will not be trivial and will probably re- 
quire simulating multiple merger events which combine to form 
the halo. This may be obtained from the results of detailed cosmo- 
logical simulations. 

• Including an extinction model and possibly adding a thick disc 
to the Galaxy model. Within the Gaia project 3D extinction models 



have been deve loped which can be added to our simulations (see 
|Primmell2003) . 

• Using more sophisticated models for the dwarf galaxies. The 
TV-body simulations should include dynamical friction and possi- 
bly more particles. The simulation of the stellar properties of the 
satellite particles should be updated with a more realistic mass 
function and possibly with the inclusion of more than one stellar 
population in the dwarf galaxy. The need for the latter is supported 
by the recent results concerning the stellar population of the Sculp- 
tor dwarf g alaxy in which two distinct ancient stellar components 
were found l Tols toy etall2004 . 
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APPENDIX A: TRANSFORMING PHASE-SPACE 
COORDINATES TO GAIA DATA 

Here we describe in detail how the phase-space coordinates (r, v) 
of our model Galaxy were converted to the data in the simulated 
Gaia catalogue; the five astrometric parameters I, b, w, [it,, and fib 
and the radial velocity u ra d- The barycentric position r b = r — r© 
and velocity v b = v — vq are related to the astrometric and radial 
velocity data as follows (see e.g. lESAl^Tl Sections 1.2 and 1.5): 

r b = Ap$/w (Al) 

and: 

v b = pfi e ,A v /uj + qp b A v /w + fw rad , (A2) 

where [p q f] is the normal triad defined below, and A v — 
1000 mas pc and A v — 4.74047... yr km s~ x designate the astro- 
nomical unit in the appropriate form. Note that if the distances are 
in units of kpc, the units of parallax and proper motion are micro- 
arcsecond (/xas) and micro- arc second per year (pias yr ), 

In the galactic coordinate system the components of the nor- 
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Table Al. Radial velocity precision (in km s 1 ) as a function of stellar 
spectral type and magnitude V. These are the numbers used in our simula- 
tions, taken from Fig. 8.14 in lESAll200Cl) . 
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mal triad [p q f ] are given by: 



Px q x 
Py 1v 
Pz q z 

— sin£ 
cos£ 




sin b cos £ cos b cos £ 
- sin b sin I cos b sin I 
cos b sin b 



(A3) 



where f is the unit vector that specifies the (instantaneous) barycen- 
tric coordinate direction (£, b) 3 , and p and q are the unit vectors in 
the directions of increasing £ and b at f . The galactic coordinate 
components of r b and v b may thus be written: 




and: 











(J! 













fii*A v /zu 
fibA v /zu 

Wrad 



(A4) 



(A5) 



These equations can be inverted to calculate the astrometric 
and radial velocity data from the barycentric positions and veloci- 
ties, yielding: 



and 




A P /\r" 



R' 



(A6) 



(A7) 



where — fit*A v /zu and vt = yn,Av/vj are the velocity com- 
ponents perpendicular to the line of sight to the star, and R' indi- 
cates the transpose of R. 

Once the astrometric data have been calculated from r b and 
v b they are transformed to the Ecliptic coordinate system before 
adding the astrometric errors. This is done to simulate the fact that 
the astrometric errors will vary over the sky as a function of the 
Ecliptic latitude (3. As a consequence, upon rotation back to the 
Galactic reference frame, correlations are introduced between the 
errors in I and b, and fig* and fib- The transformation to the Ecliptic 
coordinate system is done through a rotation matrix as explained 



J Note that the vector f and the barycentric coordinate direction are not 
strictly the same, the latter includes the effect of stellar proper motion. 
However as we are here only interested in the instantaneous positions and 
motions of the stars at the central catalogue epoch we need not concern our- 
selves with this distinction(for more details see ESA 1997, Vol. 1 Sections 
1.2 and 1.5) 



in Section 1.5.3 o f the introdu ctory volume to the Hipparcos and 
Tycho Catalogues lESAll99?t) . 

The errors in the astro metry are then added as follows (see 
|Perrvmanl2002t IeS A|2 000 ) . The error on the parallax is given by: 

a ro ~ (7 + 105z + 1.3z 2 + 6x KT 10 z 6 ) 1/2 

x[0.96 + 0.04(V-/)], (A8) 

where z = 10 0,4 ' G_15 ' and G is the broad band Gaia magnitude 
(defined by the mirror reflectivity and CCD quantum efficiency 
curves) which can be calculated from the approximate transforma- 
tion: 



G= V + 0.51 - 0.50 x vU6 + (V - / - 0.6) 2 
-0.065 X (V — I — 0.6) 2 . 



(A9) 



For the mean position and proper motion errors <jo and <r M the fol- 
lowing mean relations are used: 



(Jo = 0.87(7^ 



0.75a c 



(A10) 



Finally, the variations of the errors as a function of (3 are taken from 
Table 8.3 in lESAlfcOOol) . 

The error s on the radial velocity have been taken from fig. 8.14 
in lESAl 12000). which shows the radial velocity accuracy for a num- 
ber of values of V for both hot (OBA) and cool (FGKM) stars. The 
numbers are given in Table lAll and interpolation was used for in- 
termediate values of V. For OBA stars fainter than V = 17 and 
FGKM fainter than V — 18 no radial velocities are generated. 

Finally, the predicted observational errors for Gaia as quoted 
here reflect the estimates at the time the Gaia Concept and Tech- 
nology Study Report was published (ESA 2000). The values will 
evolve along with changes in the details of the mission design and 
implementation. 



APPENDIX B: ANGULAR MOMENTUM FROM 
ASTROMETRIC AND RADIAL VELOCITY DATA 

Because the Sun is not located at the centre of the Galaxy the re- 
lation between the angular momentum of a star and its measured 
position, parallax, proper motion and radial velocity is complicated 
by the necessity to transform from the barycentric to the Galac- 
tocentric frame. The angular momentum vector L is calculated as 
r x v = (r b + tq) x (v b + vq), which when worked out using 
Eas. <A4t and <A5t yields: 

L = L g + L® , 



where L g is given by: 

UJ 2 

and L® is given by: 



/ fib sin £ — fit* sin b cos A 

—jib cos £ — fie, sin b sin I 
\ fit, cos b J 



(Bl) 



(B2) 



/ 



V 



-Re 



A v cos b 



+ f r ad sin b 



j-) a A v cos b cos £ T r 
«©A H Vq 



(B3) 



where the term A is given by: 

A v fit* cos £ 



A =- 



A v fib sin b sin I 



+ v la d cos 6 sin £ + Vq . 



(B4) 
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L g is the 'Galactocentric' angular momentum which has no phys- 
ical meaning but corresponds to the angular momentum calculated 
by an observer at the Galactic centre. The vector L indicates the 
components of the observed angular momentum that arise from the 
transformation from a barycentric to Galactocentric frame of refer- 
ence. The norm of L s correctly evaluates to: 

|L g l = ^rV^* + ^' (B5) 

i.e., the distance to the Galactic centre multiplied by the total tan- 
gential velocity. The equation above shows how the radial velocity 
enters the angular momentum measurements and can be used to 
further study the way the astrometric errors propagate in the mea- 
sured angular momentum. 

The equations for the energy of a star can be worked out sim- 
ilarly but they depend on the model for the potential energy and 
hence we do not provide them here. 



